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Abstract. Motivated by the fact that empirical time series 
of earthquakes exhibit long-range correlations in space and 
time and the Gutenberg-Richter distribution of magnitudes, 
we propose a simple fault model that can account for these 
types of scale-invariance. It is an avalanching process that 
displays power-laws in the event sizes, in the epicenter dis- 
tances as well as in the waiting-time distributions, and also 
aftershock rates obeying a generalized Omori law. We thus 
confirm that there is a relation between temporal and spatial 
clustering of the activity in this kind of models. The fluctu- 
ating boundaries of possible slipping areas show that the size 
of the largest possible earthquake is not always maximal, and 
the average correlation length is a fraction of the system size. 
This suggests that there is a concrete alternative to the ex- 
treme interpretation of self-organized criticality as a process 
in which every small event can cascade to an arbitrary large 
one: the new picture includes fluctuating domains of coher- 
ent stress field as part of the global self-organization. More- 
over, this picture can be more easily compared with other 
scenarios discussing fluctuating correlations lengths in seis- 
mic ity. 



1 Introduction 

At the moment there is not a comprehensive explana- 
tion of the mechanisms giving rise to the complex phe- 
nomenology of earthquakes. The magnitude of each 
eart hquake is characterized by the Gutenberg-Richter (GR) 
law (IGutenberg and Richten. 119441) . which is in fact a scale- 
invariant distribution of energy release. Earthquakes are 
also long-range correlated with each other. It is indeed 



known that events are clustered in space and time dTurcotte , 
[1997; Scholz, 2002) and take place in complex fault pat- 
terns |Bonnereral|_ 2001). The Omori law of aftershocks 



rate (lUtsu et al.Lll995h is an example of the temporal cluster- 



law. The phenomenology of the distance between subse- 
quent e picenters is also characterized by power-law distri- 
butions ( Davidsen and Paczuski , 2005 ; Corral , 2006h . More- 
over, the values of magnitudes, waiting times an d locations 



of earthquakes are part of a singl e scaling picture (iBak et al. 
I2OO2I; 'Corral', '2003* '2004', '2005"). Other e xamples are given 
by Mega et al. (2003) and Davids en et al. ]|2006). bmce seis- 
micity is one of the most outstanding examples of a class of 
phenomena involving a wide range of energetic, spatial, and 
and temporal scales, it is expected that its modeling is prob- 
lematic. 

It is possible to build models based upon the phenomenol- 
ogy of earthquakes. For example, aftershock-sequence 
models req uire an assumed law of off-spring generation 



per event 



Turcotte et al.L 



JOgata . 



2007 



19881 'Helmstetter an d Sornette . 
; iLippiello et al., 2007i). These models 
can yield realistic time-series, but by construction they use 
rather than explain laws like the GR one. 

The scale-invariant distribution of earthquake sizes is re- 
produced by processes based on avalanches of stress redis- 
tribution, following the idea that there is self-organized crit- 
icaUty (SOC) (Bak, 1996; Sornette, 2000). The precursor of 
thi s concept in geophysic s has b een the slider-block model 
by iBurridge and Knopofn (119671) . It is evident from many 



2002t 



models that the mechanism of avalanches of relaxations ro- 
bustly leads to size-frequency power-laws. This behavior 
emerges from the collective organization of units that co- 
operate with very nonlinear rules, redistributing stress and 
typically dissipating it from open boundaries. 

However, it has become also clear during the last years 
that the simplest SOC models cannot reproduce other im- 
portant features of critical phenomena, usually involv- 
ing correlations between events. Models incorporat- 



ing of earthquakes, with a decay given by a scale-invariant 



ing correlated events (Olami et al.. 
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nority within the literature on SOC. These few scattered re- 
sults unfortunately have not constituted a large enough body 
for appropriately raising the issue of temporal organization 
to the attention of the scientific community. 

In this paper we show that earthquakes phenomenology 
can guide us to build self-organized models with the appro- 
priate features. In particular, we stress the importance of 
clustering events in space and time, an aspect leading us to 
develop a fault model that displays a full spectrum of power- 
law statistics (GR law, Omori law, waiting times and epicen- 
ter distances with broad distributions), not observed in pre- 
vious models. Hence, the very basic idea of SOC is in fact 
achievable. In particular, the process self-organizes the epi- 
center locations, clustering them rather than spreading them 
randomly in space, as it is frequently imposed in other simple 
models. 

A novel feature distinguishing the model we propose from 
previous ones is is the possibility to infer maximal areas of 
events from its configuration. It turns out that this model does 
not conform to the common picture associated with S OC in 
geophysics (iNature debate. 19991: lOeller et al.i Il997h . The 
idea is that every tremor can in principle cascade in a large 
event, depending on minor details of the stress field. It is 
possible that the paradigm of sandpiles has been much influ- 
ential in the consolidation of this view. Up to date, this in- 
terpretation has been a speculation, without any quantitative 
assessment of its validity. Below we show that we instead 
observe a mean correlation length limited to a given fraction 
of the whole fault, and a rich dynamical regime leading to 
complex patterns of possible slipping areas. The domains 
where avalanches can occur are not always maximal. There- 
fore, it is clear that in this model it is not possible to have 
a large earthquake at all times. We will come back to this 
point in the Section "Discussion". The next section contains 
the description of the model, while the numerical results are 
shown in section [3] 

2 Model 

The following model describes a one-dimensional fault with 
L units and with periodic boundary conditions. Each unit i 
represents the displacement hi of a plate with respect to a 
second one. Plates are sliding with respect to each other and 
thus the displacement hi corresponds to a slip accumulated 
with time. An external field ai characterizes the speed of the 
strain accumulation in the unit: At each time step a unit i, 
chosen with probability Pi exp(/3tTi), shps: 



1 



(1) 



If hi forms a high gradient with one of its neighbors j, in 
our case hi — hj > 4, a local elastic instability occurs. This 
is relaxed by allowing the two nearest-neighbor units to get 
closer. 



If this process leads to the formation of new unstable cou- 
ples they are listed and processed into a random order 
until the list is empty, filling at the same time another list 
with eventual new unstable pairs. The new list is then pro- 
cessed, and so on. The iteration of this rule leads to a final 
state in which all bonds between units are stable again. The 
whole avalanche of relaxations represents an earthquake and 
is characterized by its size (the number of single relaxations, 
corresponding to the seismic moment), by its slipping area 
(the number of sites involved at least once), and by its epi- 
center (the unit where the avalanche started). It takes place 
by definition in one time step. The waiting time between 
avalanches is then measured by the number of time steps sep- 
arating them. 

The aim of the field cr,; is to reproduce some "external" 
tectonic loading, which should be originated by the crust por- 
tions that meet at the fault. Somewhat a replaces the load- 
ing calculated explicitly with the laws of elasticity in other 



models (see for example dBen-Zion , 119961: iBen-Zion et al. 



20031)). Since earthquakes play the main role in reshaping 
the stress field in the crust, we let each ai evolve with a rule 
that couples it with the activity in the system: Every time that 
a redistribution ^ occurs, the two corresponding fields are 
set equal to their average aij = {ai+aj)/2 plus a noise term 
6 drawn at random (for each site) from the interval [— 1 , 1] Q 



and 



(3) 



The evolution of the system is thus stochastic in many as- 
pects. At the level of single redistributions involving (|2| and 
(O, one has an update of a's with random (5's. At the step ([TJ 
of forcing the system, the choice of i according to a proba- 
bility Pi is also stochastic. One can interpret the set of (Ti as 
an array of local rates. Indeed, a micro-slip (hi hi + 1) 
takes place with a rate proportional to exp(/3(Ti). 

A non-trivial regime emerges as long as /? is sufficiently 
large to lead to a persistence of the earthquake activity in ar- 
eas of the system. For (3 oo one finds a choice of the posi- 
tion to apply ([T]) that corresponds to the site with the largest 
a. This resembles an extremal dynamics for the field ai. We 
rather chose /? large but finite, such that many parts of the 
fault are likely to be active at the same time (if they are share 
similar values of ai). The evolution of the ai guarantee a 
migration of active areas as well. 

Despite the stochastic character of some of the micro- 
scopic updates, a rich phenomenology arises, with scale-free 
avalanches and with realistic interoccurrence statistics. 



3 Results 

We show results obtained by fixing /3 = 4, which is large 
enough to lead to clustering of epicenters. A preliminary 



and 



(2) 



The choice of this interval just fixes the scale of fluctuations of 
the ai's. 
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check has shown qualitatively similar results in the range 2 < 
f3 < 6. For each L, initial configurations for simplicity have 
hi = and ai = 0. To be confident that the stationary regime 
has been reached, we first run a long transient of « 10^ ^ 10^ 
time steps without collecting statistics. From time step t = 
we then collect time series composed by 2 3 x 10^ time 
steps. This constitutes a satisfactory statistics only if a large 
number of different profiles is sampled, which is the case for 
systems with < 2000 units. We can thus collect data in a 
reasonable time for systems up to this size. 

A first glance at the behavior of the model is proposed in 
Fig.d] where we plot a sample of size and location of rupture 
areas as a function of time. One can see that the activity is an 
alternation of earthquakes of several sizes, with a persistence 
in active areas. This is confirmed by a plot of the increment 
of hi with respect to the values at time t = 0: Fig.|2tb) shows 
that the increments are concentrated in the active areas. 

The statistics of several quantities turn out to be deter- 
mined by power-laws. In order to display the frequency-size 
statistics, we adopt the following definition of magnitude: 

TO = logio S 



Note that the usual prefactor 2/3 ('S cholz , 2002) in the con- 
version from seismic moment to magnitude is not suitable 
for a one-dimensional model because the area of events is 
in fact a length. In Fig. [3] one can see that the number of 
events with magnitude > to, denoted by Ny{m), seems to 
follow a GR law, Ny{m) - lO"*"™, with 6 = 1.1 ± 0.1, 
though this distribution is most likely multiscaling, as it is 
often the case in one-dimensional automata ( Kadanoff et al.L 
1989h . We postpone the exact characterization of this dis- 
tribution to future work. The distribution of slipping areas 
a instead has a clearer scaling: it develops a power-law tail 
a^"^" for increasing L, with = 1.5 (Fig|4|i, and obeys to 
standard finite-size scaling 



P(a) 



(4) 
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Fig. 1. Example of a time series for a system with L — 2048 sites: 
(a) size vs time and (b) location of rupture areas versus time. 
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Fig. 2. (a) Profiles hi corresponding to the configuration at time 
t = of Fig. [T] (black line) and at time t = 50000 (red line), and 
some intermediate stages (thin gray lines). To all curves we have 
subtracted the average h at time t = 0. (b) Difference of the same 
profiles with respect the initial one, hi{t = 0), to better visualize 
the regions where activity was concentrated in this example. 



with D = 1 and where F is a scaling function, see inset of 
Figgl 

In addition to the avalanche size and area, in this model we 
can also measure metric properties characterizing the state 
of the system between two avalanches: one is the length of 
domains of units having constant sign in the slope of hi. Each 
profile hi is indeed an alternation of domains with increasing 
h and domains of decreasing h, forming in general a non- 
trivial landscape, see Fig. 12 a). This is also a result of the 
self-organization of the process, which includes the evolu- 
tion of the ai. Also domain lengths £ have a power- law 
distribution ^ with Tg ~ 1.9, see Fig.|5] which displays 
finite-size scaling 



also with D = 1 (inset of Fig. |5]l- Since Tq < T£, there 
is more chance to observe large areas than large domains. 
On the other hand, avalanches take place within domains. 
This suggests that avalanches are repetitive and appear more 
frequently in long domains. 

Connected with the scale-invariance of domains, there is 
also a scaling of the correlation length of the stress fi — 
hi+i — hi with the system size. The correlation length can 
be read from the shape of the correlation function 



ifi+rfi) — {fi) ifi+rfi) 



inn) - {.nr 



(6) 



(5) 



where (...) means a statistical average over the sites and con- 
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Fig. 3. Gutenberg-Richter law in systems with L = 1024 and L 
2048. 



Fig. 6. Correlation function CL{r) of the stress fi for L — 256, 
512, 1024, and 2048, plotted as a function of (a) r and (b) r/L. 
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Fig. 4. Distributions of the area of avalanches, for L = 128, 256, 
512, 1024, and 2048. Their power-law tail ~ a"^'^ is highUghted 
by the dashed line. Inset: data collapse of P{a)a^°- vs a/L. 
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Fig. 5. Distribution of domain lengths £ (same L's of Fig.|4ll. The 
dashed line represents a power-law £~^'^. Inset: data collapse of 
P{l)r'^ vs l/L. Data for the shortest L = 128 are not included in 
the collapse. 




jump 

Fig. 7. Distribution of the jumps (distances between subsequent 
activities) for the same L's of Fig.|4] Their power-law tails have an 
exponent converging roughly to ~ 1 for large L. 



figuration^ It turns out that CL{r) conforms to a scaling 
function CL{r) — C{r/L), with C{. . .) independent on L, as 
shown in Fig. |6l Hence, if we define the correlation length 
as the range where Cl(?') > 0.1, we see (Fig.|6]l that it has a 
value ~ 10%L that diverges linearly with L, as one expects 
in critical systems. We will come back to this point is the 
Discussion. 

Another quantity of interest is the jump between the po- 
sition of grain addition at time t and the subsequent posi- 
tion of grain addition att + 1. The jump distributions have 
also power-law tails, with exponent converging to « —1, see 
Fig. Q This distribution is thus similar to that of distanc es 



between sub s equen t earthquakes (Davidsen and Paczuski , 
20051: ICorralL l2006h . Also the crossover to a background 
level for long jumps takes plac e at a length that is a fixed frac- 
tion the size of the catalogue ( Davidsen and Paczuski ,20()5[ 
CorralLl2006h . 



^The periodic boundary conditions imply (fi) — {fi+r) = 0. 
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Fig. 8. Distribution of waiting times, for events larger than thresh- 
olds s (L — 2048). Two power-law fits are also shown for the two 
parts of the distribution relative to s = 30000. 
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Fig. 9. Rescaled distribution of waiting times, {t^} is the mean 
waiting time between events (it depends on the threshold s). 



3 . 1 Temporal correlations 

During the last years part of the scientific debate on earth- 
quake correlations has been focu sing on the statist i cs of w ait- 
ing times between events, see aiesi a nd Maesl l2006h for 
an overview. An issue was whether SOC models can have 
avalanches correlated with each other Some models have 
waiting times between avalanches with an exponential distri- 
bution, suggesting that their events are completely uncorre- 
cted. Clearly this is an unw a nted f eatur e in models of earth- 
quakes. Recently iBak et al.l (|2002") and ' Corrall (|2003[ 12004 , 
12005) have shown that waiting times have in general a non- 
trivial scaling form in their distributions. 

In Fig. [8] we plot some waiting time distributions that we 
observe in our model, for L = 2048 and for several minimum 
thresholds s of the size. These distributions have a shape with 
a double power-law form for high thresholds, as observed in 
catalogs of regional seismicity by Corral (2003,) and in an 
aftershock-sequence model bv iLippiello et al. ( 2007b . 

In Fig.|9]there is an attempt to collapse some of these dis- 
tributions on a single curve, by rescaling the waiting times 
to scales in which their average value is 1, that is, by multi- 
plying their values by the rate of events larger than the cor- 
responding minimum thresholds. This procedu re revealed an 
interesting s caling form for real earthquakes ( CorralL 2003 



is a salient feature of seismicity, characterizin g the occur- 



'2004', '2005") (and also for solar flares, see jBaiesi et al 
.2006)): in that case one observes a nice data collapse, with 
distributions being described by a single scaling function. 
The data collapse for this model is only approximate. We 
can conclude that the power-law tails in the distributions are 
a clear indication of a non-trivial organization and clustering 
in time of the avalanches, with some missing scale-invariance 
evidenced by the thresholding procedure. 

It is also not trivial to observe aftershocks in simple mod- 
els of seismicity. Indeed, one does not always observe Omori 
decay of aftershocks in synthetic catalogs. However, this 



rence of correlated event s even for years (Utsu etal 





1995 




2005 



Zahapin et al.1 2008h . Our model does not yield time se- 
ries with patterns clearly identifiable with aftershocks se- 
quences, intended in the usual seismological sense. Never- 
theless, an Omori-like decay can be detected, confirming the 
temporal clustering evidenced by waiting time statistics. To 
visualize the Omori decay, we use a simple definition of af- 
tershocks, leaving more complicated spatio-temporal anal- 
ysis ( Shcherbato Baiesi and Paczuski, 2003, 
20051: iBaiesii I2OO6I: [Zahapinet all boOSi) for future works. 



Let us consider events with size s as main shocks (to im- 
prove the statistics, we actually consider events in a range 
[0.9s^^, l.ls*^]). Each of these events collects aftershocks 
in a time-window following its occurrence time and in- 
cluding only events of smaller size. This time window t — t^^ 
thus ends if a new event of size at least 0.9s*^ occurs. The 
averaged statistics of the rate r{t — t^^) of avalanches after 
an main event of size s*^ is shown in Fig.[TO]as a function of 
the time lag t — t^^ from the main shock, for several values 
ofs^^. 

One can see that the aftershock decays depend on s^^ and 
follow a generalized Omori decay 



rit) 



A 



(7) 



where A is a constant, t* is a characteristic time, and p is 
the exponent of the generalized decay (usually one observes 
p 1). As in real seismicity (Baiesi and Paczuski, 200^ 



20051) . the onset of the power-law decay takes place at times 
t* that increase with the size of the main event. The same 
is true for the end of the Omori decay: data in Fig. [TOlhave 
an exponential decay after the Omori regim e, as it was found 
for aftershocks ( Baiesi and Paczuski 20051) . The exponent 
p takes values ranging from « 1.3 for s*^ = 300, to « 0.5 



6 



M. Baiesi: Correlated earthquakes in a self-organized model 



-0.5 




1 1 1 1 1 


1 I 1 I 1 








-1 






^^s"^ = 10000 


? -1.5 

o 

so" 

^ -2 




s'' = 30o\ 




-2.5 








-3 




1,1,1 





3 

log,„(t-t ) 



Fig. 



10. Decay of aftershocks activity after main shocks of size 
= 300, 1000, 3000, and 10000, in a system with L = 2048. 
Dense lines are data, while dashed lines are fit according to the gen- 
eralized Omori decay Q. 



for s*^ = 10000. Its variability somewhat reflects the same 
lack of invariance for increasing thresholds manifested by 
waiting-time distributions. 

4 Discussion 

Some previous SOC models with realistic phenomenol- 
ogy are based on the m e chanism of extrem a l dy 
namics (Olami et all 1 19921; iHainzl et all 1 19991 I200C: 



Hergarten an d Neuge baueii 2002 ; ZoUer et al 



2005 



Lippiello et al.i . i2005i) . in which an earthquake starts always 
from the weakest unit. Our stochastic model shows a more 
general mechanisms giving rise to correlated events within 
SOC, which involves activity suitably clustered in space 
and time, together with scale-free redistributions of energy 
in the form of avalanches. The random aspect cannot be 
excessiv^: a load completely random in space has been 
for years the standard in several SOC cellular automata, 
maybe because it is the simplest protocol. In the field of 
seismicity this choice is not supported by phenomenological 
observations, as we know that epicenters are correlated and 
clustered. When a random lo ad was imposed, avalanches 
were found to be uncorrected ( Baiesi and Maesl 20061) . We 
thus argue that a (correct) clustering in space of events can- 
not be disentangled from the temporal clustering of events, 
both aspects being part of the same global organization in 
critical systems. 

Regardless of the lack of dissipation from open bound- 
aries, our process reaches a stationary critical regime. The 
reason is that its loading is not homogeneous and the evolu- 
tion via avalanches generates the domains over which further 



In our model, the activity spreads randomly in space with low 
/3 values. In this limit, domains shrink to exponentially short regions 
and the system loses scale-free avalanches. 



large avalanches can occur In the periodic system we have 
described, the minima of the accumulated slip profile are 
places where eventually avalanches must stop. These min- 
ima are not fixed but dynamic. 

It is important to note that the dynamics of the accumu- 
lated slip profile, with domains that evolve in time, has non- 
trivial consequences. Each domain seems to represent what 
is normally observed in canonical SOC systems with open 
boundaries dBakl Il996h . the so called "sandpiles", which 
have a profile with a single slope, from the maximum at 
a closed boundary to a minimum at an open (dissipating) 
boundary. Eventually the whole process somewhat resem- 
bles a collection of smaller homogeneous SOC systems, 
whose number and position fluctuates in time. For each con- 
figuration, the maximum correlation length should be close 
to the length of the longest domain. Interestingly, this do- 
main length is not always close to its possible maximum, 
which means that the system is often in a state incompati- 
ble with an earthquake spanning the whole fault. Moreover, 
we have seen that the range of the average correlation length 
is a fraction of the system size. On the one side, this says that 
we have to reconsider the typical value of correlation ranges 
upon change of scale of the whole system. Provided that we 
can meaningfully isolate an area from the rest of the crust, 
on the other hand, we can expect a finite mean correlation 
length within it. 

Hence, our model does not reproduce a popular picture 
associated with SOC, invoking a continuous state of "max- 
imal" criticalit y in the crust due to a n eventual infinite cor- 
relation length ([Nature debate. 1999h . According to this pic- 
ture, earthquakes are inherently unpredictable in size, space 
and time because their cascade to large events depends on 
minor details of the stress field. This point has been used, 
for example, by G eller et al.. (1997) to infer that earthquakes 
cannot be predicted. The validity of their argument can be 
limited by the lack of discussion about non-minor details. 
These major details in our models are those that are macro- 
scopically visible when looking at the profile of the slip field 
hi, namely the different domains. Unfortunately patterns like 
thes e are not accessible i n real measurements. Bak pointed 
out (INature debate. 1999b that an earthquake does not "know 
how large it will become". This is not incompatible with 
our point that an earthquake "knows how large it cannot be- 
come". Perhaps both aspects sho uld be taken into acc ount in 
studies on earthquake prediction dKeilis-Borok , 2002). 

Therefore, according to our results, the following scenario 
is possible: The process of self-organization in seismicity, 
due to the slow load of the crust and its fast relaxation via 
earthquakes, converges to a dynamical SOC regime, with rise 
and fall of patterns of strongly correlated stress. These pat- 
terns may be associated with (local) fluctuating correlation 
lengths. 

One could also have coexistence of SOC and other mecha- 



nisms (ISammis and Sornettell2002h . A previous SO C model 
with a heterogeneous fixed pattern of faults (.Huang et al 
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19981) has a behavior consistent with the hypothesis that 
the approach to large earthquakes is described by a critical- 
point picture dJaume and Svkes , Il999l: l ISammis and Sornettel 



20021) . with a finite-time singularity of Be nioff strain release 



and a divergence of a corre lation length dZoUer and Hainzl , 



2002t IZaliapin et al. . I2OO2I ). We have not investigated this 



point in our model yet, though it seems that its dynamics 
does not break all the correlations after a large earthquake. 
Indeed, a large slip along a domain lowers the total energy 
stored in the system, and eventually shifts the domain range 
of some units, but the domain itself should be ready for 
similar earthquakes without too much effort. However, an 
eventual merging with other coherent domains might lead 
to an increase of the correlation length i n the area, with a 
possible connection with previous studies ('jaume and Svkes", 
1999: Sammis and Sornette. 2002; ZoUer and Hainzl. 2002,; 



Zaliapin et al.L l2002h . In any case, the stationary regime 



of our model appears to be di fferent fro m that of intermit- 
tent c riticality (IBen-Zion et al.l 12003; Bowman and Sammis , 



20041) . in which every large event drives the system far from 



criticality, which is then slowly restored by the dynamics. 



5 Conclusions 

We have shown that it is possible to build stochastic pro- 
cesses with self-organized criticality that reproduce several 
power-laws found in earthquake statistics, like the GR law, 
the generalized Omori laws, the waiting-time distributions, 
and the distributions of distances between subsequent events. 
The robust scale-invariant statistics generated by avalanches 
is the leading principle of the study. We have stressed that 
it is important that the process generates activity clustered in 
space for eventually obtaining a clustering in time of events. 
In our model, contrary to previous examples, the clustering 
takes place even if the system is stochastic, showing that 
a moderate degree of randomness can be tolerated in SOC 
models with spatio-temporal correlation. 

Our findings are contrary to a constant complete unpre- 
dictability of event sizes, even if SOC is one of the main 
mechanisms acting to generate the complexity of seismic- 
ity. The point is that every SOC system can have a finite 
size. We have described a system displaying domains that 
resemble a fluctuating collection of canonical SOC cellular 
automata. Interestingly, domains limit each other and their 
boundaries constitute the points where avalanches eventually 
must stop. The size of each domain quantifies locally the cor- 
relation length, which is thus a quantity fluctuating in space 
and time. As a result, the mean correlation length diverges 
with the system size, as expected, but it occupies only a fi- 
nite fraction of the system. We have been able to visualize 
these features thanks to the simplicity of the one-dimensional 
model. 

The oversimplified process that we have discussed is in- 
spired by the phenomenology of earthquakes and tries to en- 



capsulate it, but clearly it is a geophysical model in an em- 
bryonic stage. Hopefully the results and discussion we have 
presented provide new ideas that will be useful for build- 
ing models grounded on laws of geophysics and elasticity 
of solids, which still preserve the ability to reproduce earth- 
quakes phenomenology. With models of this kind, for ex- 
ample, it would be interesting to see if creeping sections of 
faults can play the role of domain boundaries in the sense 
discussed in this paper. 
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